Figure 05 Sulfonylurea (SFU)

knitr::opts_chunk$set(
  warning = FALSE,
  error = TRUE, 
  echo = TRUE, 
  fig.width = 10, 
  fig.height = 7)
library("magrittr")
library("data.table")
library("ggplot2"); packageVersion("ggplot2")
## [1] '3.3.5'
library("patchwork"); packageVersion("patchwork")
## [1] '1.1.1'
library("ggbeeswarm"); packageVersion("ggbeeswarm")
## [1] '0.6.0'
library("pspline"); packageVersion("pspline")
## [1] '1.0.18'
theme_set(
  theme_bw() +
    theme(
      panel.grid = element_blank(), 
      axis.ticks = element_line(size = 0.25),
      strip.background = element_blank(),
      strip.text.y = element_text(angle = 0)
    )
)
scaleColorManualTreatment <-
  scale_color_manual(
    values = 
      c(
        placebo = "gray50",
        wbf10 = "darkblue",
        wbf11 = "darkgreen"
      )
  )
scaleFillManualTreatment <-
  scale_fill_manual(
    values = 
      c(
        placebo = "gray50",
        wbf10 = "darkblue",
        wbf11 = "darkgreen"
      )
  )
studyArmPretty <- 
  c(placebo = "Placebo", 
    wbf10 = "WBF-010", 
    wbf11 = "WBF-011")
scaleColorTreatmentManPretty <-
  scale_color_manual(
    values = 
      c(
        "Placebo" = "gray50",
        "WBF-010" = "darkblue",
        "WBF-011" = "darkgreen"
      )
  )
scaleFillTreatmentManPretty <-
  scale_fill_manual(
    values = 
      c(
        "Placebo" = "gray50",
        "WBF-010" = "darkblue",
        "WBF-011" = "darkgreen"
      )
  )
scfaNamePretty <-
  c("Acetic acid" = "Acetate", 
    "Propionic acid" = "Propionate", 
    "Butyric acid" = "Butyrate")
vecOrderScfa <- 
  c("Acetic acid", "Propionic acid", "Butyric acid") %>% 
  rev()

Load data

tabSfu <- readRDS(params$tabSfu)
tabGluCtl <- readRDS(params$tabGluCtl)
# Per protocol
tabPp <- readRDS(params$tabPp)
tabPpOrigFormat <- readRDS(params$tabPpOrigFormat)

join

tabSfuSimple <-
  tabSfu %>% copy %>% 
  .[, 
    .(
      SulfonylureaDetected = any(SulfonylureaDetected),
      SUL_TREAT = any(SUL_TREAT=="Yes"), 
      SUL_CONT = any(SUL_CONT=="Yes")
    ), 
    by = "Subject"]
tabSfuGluCtl <- 
  tabSfuSimple[tabGluCtl, on = "Subject"] %>% copy()

SFU Detection Accuracy

How well does untargeted plasma data recall (in the data sense, sensitivity) the SFU usage of subjects?

# Confusion matrix between the two clinical records for SFU drug use
tabSfuGluCtl[, .(SUL_TREAT, SUL_CONT)] %>% table()
##          SUL_CONT
## SUL_TREAT FALSE TRUE
##     FALSE    44    0
##     TRUE      3   10
# So SUL_TREAT is a very mild superset, probably of the form
# 'has this subject ever been treated with SFU',
# while `SUL_CONT` is, 'is this subject currently taking SFU continuously'
tabSfuGluCtl[, .(SulfonylureaDetected, SUL_TREAT)] %>% table() 
##                     SUL_TREAT
## SulfonylureaDetected FALSE TRUE
##                FALSE    36    2
##                TRUE      8   11

Note Subject 48 was not included in plasma measurements.

tabSfu[(Subject == "48")]
## Empty data.table (0 rows and 10 cols): treatment,SUL_TREAT,SUL_CONT,Subject,Event,Tricode...
tabSfu$Subject %>% uniqueN()
## [1] 57
tabPpOrigFormat[(Subject == "48"), .(Subject, SUL_TREAT, SUL_CONT)]
##    Subject SUL_TREAT SUL_CONT
## 1:      48        No       No
tabSfuGluCtl[(Subject == "48"), c("SUL_TREAT", "SUL_CONT") := FALSE]
# Add a variable indicating whether SFU was reasonably 'expected' to be found...
tabSfuGluCtl[, Expected := (SUL_TREAT | SUL_CONT)]
tabSfuGluCtl[, DetectedOrExpected := (Expected | SulfonylureaDetected)]

Expected v. Detected (in plasma)

tabSfuGluCtl[, .(Expected, SulfonylureaDetected)] %>% table()
##         SulfonylureaDetected
## Expected FALSE TRUE
##    FALSE    36    8
##    TRUE      2   11
tabSfuGluCtl[, .(treatment, Expected, SulfonylureaDetected)] %>% table()
## , , SulfonylureaDetected = FALSE
## 
##          Expected
## treatment FALSE TRUE
##   placebo    11    0
##   wbf10      11    2
##   wbf11      14    0
## 
## , , SulfonylureaDetected = TRUE
## 
##          Expected
## treatment FALSE TRUE
##   placebo     1    4
##   wbf10       3    4
##   wbf11       4    3
# Number of subjects per detected SFU molecule
tabSfu %>% 
  melt.data.table(
    id.vars = c("treatment", "Subject", "Event"), 
    measure.vars = c("Glimepiride", "glipizide", "glyburide"), 
    variable.name = "SFU", 
    value.name = "Signal") %>% 
  .[!is.na(Signal)] %>% 
  .[, .(N = uniqueN(Subject)), by = "SFU"]
##            SFU  N
## 1: Glimepiride 10
## 2:   glipizide  9
## 3:   glyburide  5
# Summarize agreement between timepoints
tabSummarizeAgreement <-
  tabSfu %>% 
  .[, .(
    anySfuDetected = sum(SulfonylureaDetected),
    sfuResultAgrees = SulfonylureaDetected %>% equals(SulfonylureaDetected[1]) %>% all(),
    glimepirideAgrees = is.na(Glimepiride) %>% equals(is.na(Glimepiride[1])) %>% all(),
    glipizideAgrees = is.na(glipizide) %>% equals(is.na(glipizide[1])) %>% all(),
    glyburideAgrees = is.na(glyburide) %>% equals(is.na(glyburide[1])) %>% all(),
    sfusDetected = 
      c(ifelse(is.na(Glimepiride), NA_character_, "Glimepiride"),
        ifelse(is.na(glipizide), NA_character_, "glipizide"),
        ifelse(is.na(glyburide), NA_character_, "glyburide")) %>% 
      unique() %>% sort() %>% 
      paste0(collapse = ", ")
  ), 
  by = "Subject"]
# For every subject, both their blood samples agrees on SFU+ or SFU-
tabSummarizeAgreement[, .(sfuResultAgrees)] %>% table()
## .
## TRUE 
##   57
# Zoom-in on the examples where Baseline & Week12 
# disagree on which SFU drug
tabSfu[tabSummarizeAgreement[grep(",", fixed = TRUE, x = sfusDetected)]$Subject, 
       on = "Subject"] %>% 
  .[, .(treatment, Subject, Event, SUL_TREAT, 
        Glimepiride, glipizide, glyburide)] %>% 
  knitr::kable()
treatment Subject Event SUL_TREAT Glimepiride glipizide glyburide
placebo 68 Baseline Yes NA 3902473 NA
placebo 68 Week12 Yes 1491146 6712735 NA
wbf10 69 Baseline No NA NA 2053069
wbf10 69 Week12 No 5090322 NA NA
wbf10 119 Baseline Yes 573746 NA 48219
wbf10 119 Week12 Yes 6463266 NA 76898
wbf11 102 Baseline Yes 55790 NA NA
wbf11 102 Week12 Yes NA 4754458 NA
wbf11 75 Baseline Yes NA NA 343940
wbf11 75 Week12 Yes 2313739 NA NA

So all subject timepoints agree that a subject has at least one SFU Only 5 of 19 subjects with detected SFU have any disagreement or multiple SFU detected. For Subject 68 the ‘disagreement’ is trivial, in that glipizide was detected in both, but glimepiride was also detected at Week 12

tabSfu[(SulfonylureaDetected), uniqueN(Subject)]
## [1] 19

Fig. 4a: SFU and HbA1c

pSfuA1c <- NULL
vecPaletteSfu <- 
  c(`TRUE` = "orange", 
    `FALSE` = "black", 
    `NA` = "gray")
tabSfuGluCtlEmbedLegend <-
  tabSfuGluCtl %>% copy %>% 
  # only write the labels on the placebo group (left)
  .[(treatment == "placebo")] %>% 
  .[, Treatment := studyArmPretty[treatment]] %>% 
  .[, thisLegend := NA_character_] %>% 
  .[(deltaA1c == max(deltaA1c[which(DetectedOrExpected == TRUE)], na.rm = TRUE)), 
    thisLegend := "SFU\nused"] %>% 
  .[(deltaA1c == min(deltaA1c[which(DetectedOrExpected == FALSE)], na.rm = TRUE)), 
    thisLegend := "No SFU used"] %>% 
  .[!is.na(thisLegend)]
pSfuA1c <-
  tabSfuGluCtl %>% copy %>% 
  .[, Treatment := studyArmPretty[treatment]] %>% 
  .[, dummyForTitle := "\u0394 HbA1c [%]"] %>% 
  ggplot(
    mapping = aes(
      x = Treatment, 
      y = deltaA1c
    )
  ) +
  facet_wrap(~dummyForTitle) +
  geom_hline(yintercept = 0.0, size = 0.25, color = "darkgray") +
  stat_summary(
    width = 0.65,
    alpha = 1, 
    size = 0.2,
    mapping = aes(color = DetectedOrExpected), 
    geom = "crossbar",
    fun = mean, 
    fun.min = mean, 
    fun.max = mean,
    data = tabSfuGluCtl %>% copy %>% 
      .[, Treatment := studyArmPretty[treatment]] %>% 
      # Don't crossbar the lone unmeasured subject
      .[!is.na(DetectedOrExpected)]
  ) +
  geom_beeswarm(
    shape = 21,
    mapping = aes(color = DetectedOrExpected, fill = DetectedOrExpected),
    groupOnX = TRUE,
    cex = 4,
    size = 2, 
    stroke = 0.1,
    alpha = 1
  ) +
  scale_fill_manual(
    name = "SFU:",
    values = vecPaletteSfu
  ) +
  scale_color_manual(
    name = "SFU:",
    values = vecPaletteSfu
  ) +
  # Label highest and lowest points as in-chart legend
  ggrepel::geom_label_repel(
    arrow = arrow(length = unit(0.015, "npc")),
    seed = 711,
    force = 3,
    force_pull = 0.25,
    box.padding = 2,
    point.padding = 0.25,
    nudge_x = -0.2,
    nudge_y = c(-1, 1),
    data = tabSfuGluCtlEmbedLegend,
    mapping = aes(label = thisLegend, 
                  color = DetectedOrExpected),
    min.segment.length = 0,
    size = 2.25
  ) +
  theme(
    title = element_blank(),
    strip.text = element_text(size = 14, hjust = 0),
    strip.background = element_blank(),
    axis.text.y = element_text(size = 14),
    axis.text.x = element_text(size = 16),
    axis.title.y = element_blank(),
    axis.title.x = element_blank(),
    axis.ticks.y = element_line(size = 0.25),
    axis.ticks.x = element_blank(),
    legend.position = "none",
    panel.grid = element_blank()
  )
pSfuA1c

  • Which SFU was used by that single +SFU A1c ‘responder’? … Glipizide.
  • Are they the only Glipizide user? … No.  Looks like five of seven SFU-using subjects in WBF-011 were using Glipizide. So, this is only partially informative, in that it is the SFU to which the formulation strains were least sensitive during the in vitro monoculture assay.
  • The delta-A1c ordered table below does not look stratified by SFU.
tabSfuGluCtl %>% 
  .[(DetectedOrExpected & treatment == "wbf11")] %>% 
  .[, .(Subject, deltaA1c)] %>% 
  unique() %>% 
  tabSfu[., on = "Subject", nomatch = 0] %>% 
  .[(treatment == "wbf11")] %>% 
  .[, .(Subject, Event, Glimepiride, glipizide, glyburide, deltaA1c)] %>% 
  # unique() %>% 
  setorder(deltaA1c) %>% 
  knitr::kable()
Subject Event Glimepiride glipizide glyburide deltaA1c
38 Baseline NA NA NA -0.5
38 Week12 NA 6144736.00 NA -0.5
38 Week12 NA 4480478.00 NA -0.5
117 Baseline NA 20772356.00 NA 0.2
117 Week12 NA 1348684.00 NA 0.2
46 Baseline NA 601454.00 NA 0.4
46 Week12 NA 25086582.00 NA 0.4
102 Baseline 55790 NA NA 0.8
102 Week12 NA 4754458.00 NA 0.8
75 Baseline NA NA 343940 0.8
75 Week12 2313739 NA NA 0.8
18 MP07 NA 93672.21 NA 0.8
18 MP18 NA NA NA 0.8
67 Baseline NA NA NA 1.0
67 Week12 624306 NA NA 1.0

Fig. 4b: Re-estimated glucose control endpoints

tabSfuGluCtl %>% 
  .[!(DetectedOrExpected)] %>% 
  .[, mean(deltaA1c, na.rm = TRUE), by = "treatment"] %>% 
  .[, V1[treatment == "wbf11"] - V1[treatment == "placebo"]]
## [1] -0.9104895
# Original
lm(formula = deltaA1c ~ 0 + treatment, data = tabSfuGluCtl)
## 
## Call:
## lm(formula = deltaA1c ~ 0 + treatment, data = tabSfuGluCtl)
## 
## Coefficients:
## treatmentplacebo    treatmentwbf10    treatmentwbf11  
##           0.3938            0.1700           -0.2100
# Clinical record only
lm(formula = deltaA1c ~ 0 + treatment, data = tabSfuGluCtl[!(Expected)])
## 
## Call:
## lm(formula = deltaA1c ~ 0 + treatment, data = tabSfuGluCtl[!(Expected)])
## 
## Coefficients:
## treatmentplacebo    treatmentwbf10    treatmentwbf11  
##           0.3833            0.1429           -0.3529
# Clinical record or direct detection
lm(formula = deltaA1c ~ 0 + treatment, data = tabSfuGluCtl[!(DetectedOrExpected)])
## 
## Call:
## lm(formula = deltaA1c ~ 0 + treatment, data = tabSfuGluCtl[!(DetectedOrExpected)])
## 
## Coefficients:
## treatmentplacebo    treatmentwbf10    treatmentwbf11  
##           0.3182            0.1500           -0.5923
lm(formula = deltaA1c ~ DetectedOrExpected + treatment, data = tabSfuGluCtl)
## 
## Call:
## lm(formula = deltaA1c ~ DetectedOrExpected + treatment, data = tabSfuGluCtl)
## 
## Coefficients:
##            (Intercept)  DetectedOrExpectedTRUE          treatmentwbf10  
##                 0.2338                  0.5118                 -0.2552  
##         treatmentwbf11  
##                -0.6229
extract_ttest_wingrp = function(ttestOut){
  data.table(
    pvalue = ttestOut$p.value[1],
    statistic = ttestOut$statistic[1],
    estimate = ttestOut$estimate[1],
    confIntLo = ttestOut$conf.int[1],
    confIntHi = ttestOut$conf.int[2]
  )
}
extract_ttest_bwgrp = function(ttestOut){
  data.table(
    pvalue = ttestOut$p.value[1],
    statistic = ttestOut$statistic[1],
    estimate = ttestOut$estimate[1] - ttestOut$estimate[2],
    confIntLo = ttestOut$conf.int[1],
    confIntHi = ttestOut$conf.int[2]
  )
}
tabDataForTests <- tabTests <- NULL
tabDataForTests <- 
  rbindlist(
    use.names = TRUE, 
    fill = TRUE, 
    idcol = "inclusionCriteria",
    l = list(
      perProtocol = tabSfuGluCtl %>% copy,
      noSfuClinicOnly = tabSfuGluCtl[(Expected == FALSE)] %>% copy,
      noSfuClinicOrPlasma = tabSfuGluCtl[(DetectedOrExpected == FALSE)] %>% copy
    )
  ) %>% 
  # Select the variables to carry forward into tests
  .[, .(deltaA1c, deltaAucTot, deltaAucInc, 
        treatment, Subject, inclusionCriteria)] %>% 
  # pivot longer so that can loop once
  melt.data.table(
    id.vars = c("treatment", "Subject", "inclusionCriteria"), 
    measure.vars = c("deltaA1c", "deltaAucTot", "deltaAucInc"), 
    variable.name = "gluCtlMeasure", 
    value.name = "gluCtlValue")
# Within-group t-tests
tabTestsWinGrp <- 
  tabDataForTests %>% 
  .[, .(ttest = list(t.test(x = gluCtlValue, 
                            conf.level = 0.95, 
                            alternative = "two.sided"))), 
    by = .(gluCtlMeasure, inclusionCriteria, treatment)] %>% 
  .[, try({extract_ttest_wingrp(ttest[[1]])}, silent = TRUE), 
    by = .(gluCtlMeasure, inclusionCriteria, treatment)] %>% 
  setnames("pvalue", "pWinGrp") %>% 
  setnames("estimate", "estWinGrp")
# tabTestsWinGrp %>% show()
# Between-group WBF-011, t-tests
tabTestsBwGrp <- 
  tabDataForTests %>% 
  .[, .(ttest = list(t.test(x = gluCtlValue[(treatment == "wbf11")],
                            y = gluCtlValue[(treatment == "placebo")],
                            conf.level = 0.95, 
                            alternative = "two.sided"))), 
    by = .(gluCtlMeasure, inclusionCriteria)] %>% 
  .[, try({extract_ttest_bwgrp(ttest[[1]])}, silent = TRUE), 
    by = .(gluCtlMeasure, inclusionCriteria)] %>% 
  setnames("pvalue", "pBwGrp") %>% 
  setnames("estimate", "estBwGrp")

# Join the two for plotting
tabTests <- NULL
tabTests <-
  tabTestsBwGrp %>% 
  .[, .(gluCtlMeasure, inclusionCriteria, pBwGrp, estBwGrp)] %>% 
  .[tabTestsWinGrp, on = .(gluCtlMeasure, inclusionCriteria)]
charInclusionPretty <-
  c(
    "perProtocol" = "Per\nProtocol", 
    "noSfuClinicOnly" = "No SFU: Clinic", 
    "noSfuClinicOrPlasma" = "No SFU:\nClinic or Plasma")
tabTests %>% copy %>% 
  .[, Treatment := studyArmPretty[treatment]] %>% 
  .[, Inclusion := factor(charInclusionPretty[inclusionCriteria], 
                          levels = charInclusionPretty)] %>% 
  ggplot(aes(Treatment, estWinGrp, color = Treatment, fill = Treatment)) +
  facet_grid(gluCtlMeasure ~ Inclusion, scales = "free_y") +
  geom_hline(yintercept = 0.0, size = 0.25) +
  # Enforce symmetric y-axis ranges
  geom_blank(mapping = aes(Treatment, -estWinGrp)) +
  # The within-group estimates + C.I.s
  geom_pointrange(
    alpha = 1.0,
    size = 1,
    fatten = 3,
    shape = 23,
    stroke = 0,
    # position = position_dodge(width = 0.6),
    mapping = aes(ymin = confIntLo, ymax = confIntHi)
  ) +
  scaleColorTreatmentManPretty +
  scaleFillTreatmentManPretty +
  theme(
    legend.position = "none", 
    axis.text.x = element_text(
      size = 8,
      angle = 35, 
      hjust = 1, 
      vjust = 1))

unicodeDelta = "\u0394"
gluCtlMeasurePretty = c(
  deltaA1c = paste0(unicodeDelta, " HbA1c [%]"),
  deltaAucInc = paste0(unicodeDelta, " AUC [incremental]"),
  deltaAucTot = paste0(unicodeDelta, " AUC [total]")
)
tabTestsBwGrp$gluCtlMeasure %>% unique()
## [1] deltaA1c    deltaAucTot deltaAucInc
## Levels: deltaA1c deltaAucTot deltaAucInc
# Prepare plot
tabTestsBwGrpPlot <-
  tabTestsBwGrp %>% copy %>% 
  .[, GlucoseControl := factor(gluCtlMeasurePretty[gluCtlMeasure], 
                               levels = gluCtlMeasurePretty)] %>%   
  .[, Inclusion := factor(charInclusionPretty[inclusionCriteria], 
                          levels = charInclusionPretty)] %>% 
  .[, py := Inf, by = .(GlucoseControl)]
# annotation text size
sizeTextAnnotate <- 2.15
# Define plot
pTestsBwGrp <-
  tabTestsBwGrpPlot %>% 
  ggplot(aes(Inclusion, estBwGrp)) +
  # Align origin
  geom_blank(mapping = aes(y = -0.2 * estBwGrp)) +
  coord_flip() + 
  facet_wrap(~GlucoseControl, 
             ncol = 1, 
             scales = "free_x", 
             strip.position = "top") +
  # The mean
  geom_col(fill = "gray", color = NA, width = 0.5) +
  # Highlight origin
  geom_hline(yintercept = 0.0, size = 0.2, color = "black") +
  # Emphasize the per-protocol mean
  geom_hline(
    linetype = 3,
    data = tabTestsBwGrp %>% 
      copy %>% 
      .[(inclusionCriteria == "perProtocol")] %>% 
      .[, GlucoseControl := factor(gluCtlMeasurePretty[gluCtlMeasure], 
                                   levels = gluCtlMeasurePretty)] %>%  
      .[, Inclusion := factor(charInclusionPretty[inclusionCriteria], 
                              levels = charInclusionPretty)],
    mapping = aes(yintercept = estBwGrp), 
    size = 0.25) +
  # The between-group estimates + C.I.s
  geom_pointrange(
    color = "gray50",
    fill = "gray50",
    alpha = 1,
    size = 0.3,
    fatten = 2,
    shape = 23,
    stroke = 0,
    # position = position_dodge(width = 0.6),
    mapping = aes(ymin = confIntLo, ymax = confIntHi)
  ) +
  # p-value
  geom_text(
    size = sizeTextAnnotate,
    vjust = -0.5,
    hjust = 1.1,
    mapping = aes(
      y = py,
      label = round(pBwGrp, digits = 3) %>% paste("p =", .)
    )
  ) +
  # estimate
  geom_text(
    data = tabTestsBwGrpPlot[(gluCtlMeasure == "deltaA1c")],
    size = sizeTextAnnotate,
    nudge_x = 0.3,
    nudge_y = 0.005,
    mapping = aes(
      y = estBwGrp - (abs(confIntLo - estBwGrp)/10),
      label = round(estBwGrp, digits = 1) %>% format(digits = 1, justify = "right")
    )
  ) +
  ggtitle("Between-group comparisons [WBF-011 - Placebo]") +
  theme(
    plot.title = element_text(size = 10, color = "gray40"),
    axis.title = element_blank(),
    axis.text.x = element_text(hjust = 0.5, size = 10),
    strip.text = element_text(hjust = 0, size = 14)
  )
pTestsBwGrp

Fig. 4c: SFU Sensitivity, in vitro Monoculture

tabInvitroSfuSensitivity <- NULL
tabInvitroSfuSensitivity <-
  readRDS(file = params$tabInvitroSfuSensitivity)
# Define the experimental replication group.
# The same well coordinate appears many times,
# so we need to include the other 'run' identifiers
# that indicate a physically unique growth curve
tabInvitroSfuSensitivity[, repgrp := paste0(Strain, Round, plateType, WellID)]

Convert mass concentrations to millimoles per liter (mM). Molecular weights for each molecule:

molMass <-
  c(
    Glimepiride = 490.617,
    Glipizide = 445.536,
    Glyburide = 494.004,
    Sulfadiazine = 250.278
  )

Apply to table.

tabInvitroSfuSensitivity[, soluteMoleMass := molMass[Solute]]
tabInvitroSfuSensitivity[, ConcentrationMillimolar := 
                           1000 * Concentration / soluteMoleMass]
tabInvitroSfuSensitivity[is.na(ConcentrationMillimolar),
                         ConcentrationMillimolar := 0.0]

Quality Control EDA

Check controls

tabInvitroSfuSensitivity %>% 
  .[(Solute == "None")] %>% 
  .[(SolventPercent <= 3)] %>% 
  copy %>% setorder(repgrp, Hours) %>% 
  ggplot(aes(Hours, OD, group = repgrp)) +
  facet_grid(Strain ~ Inoculation) +
  geom_path(size = 0.1) +
  scale_color_brewer(palette = "Set2") +
  theme(
    panel.border = element_blank(),
    panel.grid = element_blank(),
    legend.position = "bottom") +
  ggtitle("Controls only. Spaghetti dump")

# Compare number of replicates making it beyond a trivial OD
# (this is a crude indicator of "growth", one which we'll refine later on)
# Clearly the inoculated conditions surpass an OD threshold 
# much more frequently than the NIC.
tabInvitroSfuSensitivity %>% 
  .[(Solute == "None")] %>% 
  .[(SolventPercent <= 3)] %>% 
  copy %>% 
  .[, .(frac = uniqueN(repgrp[(OD > 0.2)]) / uniqueN(repgrp)), 
    by = .(Strain, Inoculation)] %>% 
  ggplot(aes(Inoculation, frac)) +
  facet_wrap(~Strain, nrow = 1) +
  geom_col()

Effect of DMSO.

SolventPercent is final DMSO percentage in the growth media, the carryover from preparation of the dilution series stock. (SFUs have very low solubility in water, but manageable in DMSO).

tabInvitroSfuSensitivity %>% 
  .[(Solute == "None")] %>% 
  .[(Inoculation == "Inoculated")] %>% 
  copy %>% 
  ggplot(aes(Hours, OD, group = repgrp, 
             color = SolventPercent)) +
  facet_grid(Strain ~ SolventPercent) +
  geom_path(size = 0.1) +
  theme(
    panel.border = element_blank(),
    panel.grid = element_blank(),
    legend.position = "bottom") +
  ggtitle("Inoculated controls only. Spaghetti dump")

A little hard to read, so can convert to fraction above threshold as with the inoculation comparison.

tabInvitroSfuSensitivity %>% 
  .[(Solute == "None")] %>% 
  .[(Inoculation == "Inoculated")] %>% 
  .[, .(frac = uniqueN(repgrp[(OD > 0.2)]) / uniqueN(repgrp)), 
    by = .(Strain, SolventPercent)] %>% 
  setorder(Strain, SolventPercent) %>% 
  ggplot(aes(SolventPercent, frac)) +
  facet_grid(Strain ~ .) +
  geom_path(size = 2) +
  geom_col() +
  theme(
    panel.border = element_blank(),
    panel.grid = element_blank(),
    legend.position = "bottom") +
  ggtitle("Inoculated controls only.")
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

For these strains, growth tends to fall off of a cliff above 3% DMSO.

Drop the exploratory high DMSO concentrations that are not useful for comparative controls in the analysis

tabInvitroSfuSensitivity %>% nrow()
## [1] 298593
tabInvitroSfuSensitivity <- 
  tabInvitroSfuSensitivity[(SolventPercent <= 3)]
tabInvitroSfuSensitivity %>% nrow()
## [1] 296685

Extract growth window

################################################################################
#' Mask raw growth plate OD data to just the main growth phase 
#' 
#' While other features might be interesting or useful in special cases,
#' essentially all growth-curve experiments need to estimate parameters
#' related to the exponential(ish) period of growth.
#' The goal of this function is to apply rules that do a pretty good job
#' of finding this region of the data most of the time,
#' even in the presence of some funny business in the curves.
#' 
#' @param odtab A data.table with at least the following columns:
#'  `Hours`, `OD`, `WellID`. 
#' @param chunkMinHours Minimum duration for a valid growth phase, in hours.
#' @param minRangeOd Minimum range for a valid growth phase, in units of OD.
#' 
#' @return 
#' A named list with two tables.
#' 
#' - `growth` -- a [data.table] with only entries determined 
#' to be part of the main growth phase.
#' - `annotated` -- a [data.table] with all original entries, 
#' as well as new annotation columns that were created
#' as part of the growth-phase detection process.
#' 
#' importFrom pspline sm.spline predict.sm.spline
#' 
#' @import pspline
#' 
mask_growth_od = function(odtab, chunkMinHours = 1.5, minRangeOd = 0.02){
  require(pspline)
  newtab <-
    odtab %>% copy %>%
    # Enforce order on Hours
    setorder(Hours) %>%
    # Sanitize values for smoothing-spline 
    .[is.finite(Hours)] %>% 
    .[is.finite(OD)] %>% 
    # Smoothed-OD for empirical derivative is the way to go
    # otherwise noise can give you spurious changes in sign
    .[, odSpline := predict(sm.spline(Hours, OD), Hours, 0)] %>%
    .[, diff1 := predict(sm.spline(Hours, OD), Hours, 1)] %>%
    .[, diff2 := predict(sm.spline(Hours, OD), Hours, 2)] %>%
    # Note if derv is pos
    .[, PositiveDiff1 := diff1 > 0.0] %>%
    # Label "chunks" based on if we're in a change from negative/zero to positive
    .[, previous := shift(PositiveDiff1, type = "lag", fill = TRUE)] %>%
    .[, AddChunkCounter := PositiveDiff1 == TRUE & previous == FALSE] %>%
    .[, Chunk := cumsum(AddChunkCounter)] %>%
    # Compute Chunk range-difference
    .[, ChunkRangeDiff := max(OD, na.rm = TRUE) - min(OD, na.rm = TRUE), by = "Chunk"] %>%
    # Compute duration
    .[, chunkDuration := 
        max(Hours[PositiveDiff1], na.rm = TRUE) - 
        min(Hours[PositiveDiff1], na.rm = TRUE), by = "Chunk"] %>% 
    # Compute the sum of positive-diff1 (~auc with even spacing)
    .[, ChunkSumPositive := sum(diff1[diff1 > 0.0], na.rm = TRUE), by = "Chunk"]
  # Define the masked data only
  growthTab <-
    newtab %>% copy %>% 
    # Filter on duration
    .[(chunkDuration >= chunkMinHours)] %>% 
    # Consider only big-enough range diff1
    .[(ChunkRangeDiff > minRangeOd)] %>% 
    # Select the chunk with largest integral approximation for positive diff1
    .[ChunkSumPositive == max(ChunkSumPositive, na.rm = TRUE)] %>% 
    # Include only positive diff1
    .[(PositiveDiff1)]
  # For debugging purposes later,
  # return a list with both the full annotated table,
  # as well as the masked table
  return(
    list(
      growth = growthTab, 
      annotated = newtab)
  )
}
tabGrowthMask <- 
  tabInvitroSfuSensitivity %>% 
  # split(by = c("Round", "Strain", "plateType", "WellID")) %>%
  split(by = "repgrp") %>%
  lapply(FUN = function(x){
    suppressWarnings(
      # Annotate regions of the growth curve for param interpretation
      # This creates two results tables: `$growth` and `$annotated`
      mask_growth_od(x) %>% .$growth
    )
  }) %>% 
  rbindlist(use.names = TRUE, idcol = ".check")
# Sanity check
setorder(tabGrowthMask, Strain, Round, plateType, WellID, Hours)
tabGrowthMask %>% 
  .[, .SD[.N],  by = .(Strain, Round, plateType, WellID)] %>% 
  .[, .(Strain, Round, plateType, WellID, Hours, OD, odSpline)] %>% 
  .[sort(sample(.N, 10))] %>% 
  knitr::kable()
Strain Round plateType WellID Hours OD odSpline
AMUC Round1 Plain G05 20.321111 0.180 0.1818579
AMUC Round3 Plain D09 57.321111 0.323 0.3220033
BINF Round2 Plain B02 17.987778 1.087 1.0799692
BINF Round2 Plain F08 19.654444 0.150 0.1500846
BINF Round5 Plain G03 23.654444 0.858 0.8539474
CBUT Round1 Plain D10 22.321111 0.240 0.2385263
CBUT Round2 Plain B10 7.987778 1.948 1.9479442
CBUT Round3 wBCAA E04 22.654444 0.135 0.1355250
CBUT Round4 Plain F04 18.654444 0.937 0.9361028
EHAL Round1 Plain C04 23.987778 0.256 0.2565386
# # Define each time-step. Have to order first.
setorder(tabGrowthMask, Strain, Round, plateType, WellID, Hours)
# Baseline subtract the OD. 
# Since this is already masked to the presumed growth window, 
# and up-sorted by time, the first value 
# is the beginning of the growth window 
# and should be subtracted from all (by well).
tabGrowthMask[, OdSplineBaselineSubtract := odSpline - odSpline[1], 
              by = .(Strain, Round, plateType, WellID)]
# In case it dipped below zero, don't want to 
# skew later calculations (e.g. AUC, etc).
# Set to zero
tabGrowthMask[(OdSplineBaselineSubtract < 0.0), 
              OdSplineBaselineSubtract := 0.0]

tabInvitroSfuSensitivity %>% nrow()
## [1] 296685
tabGrowthMask %>% nrow()
## [1] 78863

Annotate the titration rank order within each replication group.

tabTitrationStep <-
  tabGrowthMask %>% 
  .[, .(Strain, Round, plateType, Solute, Concentration)] %>% 
  unique() %>% 
  setorder(Strain, Round, plateType, Solute, Concentration) %>% 
  .[, Titration := seq(1, .N, 1),
    by = .(Strain, Round, plateType, Solute)]
# Add back via join
tabGrowthMask <-
  tabTitrationStep[tabGrowthMask, 
                   on = .(Strain, Round, plateType, Solute, Concentration)]
# enforce order again after join.
setorder(tabGrowthMask, Strain, Round, plateType, WellID, Hours)

Summarize replicates, SFUs, design by round.

tabGrowthMask %>% 
  .[(Solute != "None")] %>% 
  .[, .(N = uniqueN(WellID)), 
    by = .(Strain, Solute, SolventPercent, Round, plateType)] %>% 
  dcast.data.table(SolventPercent + Round + plateType + Solute ~ Strain, 
                   value.var = "N", fun.aggregate = sum) %>% 
  setnames("SolventPercent", "DMSO(%)") %>% 
  knitr::kable()
DMSO(%) Round plateType Solute AMUC BINF CBEI CBUT EHAL
1 Round2 Plain Glimepiride 15 15 15 15 15
1 Round2 Plain Glipizide 10 15 15 15 15
1 Round2 Plain Glyburide 8 15 15 15 15
1 Round2 Plain Sulfadiazine 9 13 15 15 15
1 Round4 Plain Glyburide 30 0 30 30 0
1 Round4 wBCAA Glyburide 30 0 30 30 0
2 Round3 Plain Glimepiride 15 0 13 13 0
2 Round3 Plain Glipizide 15 0 15 14 0
2 Round3 Plain Glyburide 14 0 14 15 0
2 Round3 Plain Sulfadiazine 9 0 7 8 0
2 Round3 wBCAA Glimepiride 13 0 15 13 0
2 Round3 wBCAA Glipizide 12 0 15 15 0
2 Round3 wBCAA Glyburide 3 0 15 14 0
2 Round3 wBCAA Sulfadiazine 8 0 15 9 0
2 Round5 Plain Glimepiride 0 20 0 0 20
2 Round5 Plain Glipizide 0 20 0 0 20
2 Round5 Plain Glyburide 0 20 0 0 20
3 Round1 Plain Glimepiride 15 15 0 15 1
3 Round1 Plain Glipizide 15 15 0 15 13
3 Round1 Plain Glyburide 10 12 0 15 0
3 Round1 Plain Sulfadiazine 15 15 0 15 6

Show extracted portion of growth curve. Much easier to see the design by round and DMSO%. The solutes (SFUs) and their concentrations are not mapped on, but the difference in lag and spread is already apparent at a high level.

tabGrowthMask %>% 
  .[(Inoculation == "Inoculated")] %>% 
  setorder(Strain, Round, plateType, WellID, Hours) %>%
  ggplot(aes(Hours, OdSplineBaselineSubtract, 
             color = Round,
             group = repgrp)) + 
  facet_grid(Strain ~ SolventPercent, labeller = label_both) +
  geom_path(size = 0.25) +
  scale_color_brewer(palette = "Set2") +
  theme(legend.position = "bottom") +
  theme(panel.grid = element_blank())

  • BINF and EHAL were failed runs.
  • 1% DMSO is preferred where available (Rounds 2, 4).
  • Round 2 did not have BCAA.
  • Round 4 included only Glyburide.

Remake previous plot but with the obvious cuts for zooming in. - drop 3% DMSO - drop EHAL, BINF - drop sulfadiazine. It tended to be a more severe effect. More anomalies to deal with, and does not need to be in the main figure(s).

Stagger the height so that patterns/problems between rounds are discernible

tabGrowthMask %>% copy %>% 
  .[(Inoculation == "Inoculated")] %>% 
  .[(Strain %in% c("AMUC", "CBEI", "CBUT"))] %>%
  .[(Solute != c("Sulfadiazine"))] %>% 
  .[(SolventPercent <= 2)] %>% 
  setorder(Strain, Round, plateType, WellID, Hours) %>%
  # Stagger rounds on the vertical so that can contrast better
  .[, OdBgsStagger := OdSplineBaselineSubtract + 
      as.numeric(gsub("Round", "", Round))] %>% 
  ggplot(aes(Hours, OdBgsStagger, 
             color = Round,
             group = repgrp)) + 
  facet_grid(Strain ~ SolventPercent, labeller = label_both) +
  geom_path(size = 0.25) +
  scale_color_brewer(palette = "Set2") +
  theme(legend.position = "bottom") +
  theme(panel.grid = element_blank())

Fig. 4c, Curation

Zoom in on each strain and growth round.

Plot, highlight flagged wells

Function to plot flagged wells

plot_flag = function(iStrain, iRound, iTabFlag,
                     tabGm = tabGrowthMask){
  psr <- tabsr <- NULL
  # on.exit(return(list(psr, tabsr, iTabFlag, iStrain, iRound)))
  
  tabsr <-
    tabGm %>% copy %>% 
    .[(Inoculation == "Inoculated")] %>% 
    .[(Strain == iStrain)] %>% 
    .[(Solute != c("Sulfadiazine"))] %>% 
    .[(Round == iRound)] %>% 
    .[, ConcRound := round(Concentration, digits = 2)] %>% 
    setorder(Strain, Round, plateType, WellID, Hours)
  
  # Define plot
  psr <-
    tabsr %>% 
    ggplot(aes(Hours, OdSplineBaselineSubtract, 
               color = plateType,
               # alpha = Concentration,
               group = repgrp)) + 
    facet_grid(
      rows = vars(Titration),
      cols = vars(Solute)
    ) +
    geom_path(alpha = 0.25, size = 0.5) +
    geom_path(
      alpha = 1.0, 
      data = tabsr %>% 
        .[iTabFlag, on = c("Strain", "Round", "plateType", "WellID")]) +
    geom_text(
      size = 2,
      nudge_x = 2,
      mapping = aes(label = WellID), 
      data = tabsr %>% 
        .[iTabFlag, on = c("Strain", "Round", "plateType", "WellID")] %>%
        .[, .SD[which.max(Hours)], by = c("plateType", "WellID")]
    ) +
    scale_color_brewer(palette = "Set2") +
    theme(legend.position = "bottom") +
    theme(panel.grid = element_blank()) +
    ggtitle(
      label = paste(iStrain, 
                    "+ SFUs, Round", 
                    gsub("Round", "", iRound)))
  return(psr)
}

AMUC curation by round

AMUC, Round 4

tabFlagAmucR4 <-
"
Strain Round plateType WellID
AMUC Round4 wBCAA A08
AMUC Round4 wBCAA A12
AMUC Round4 Plain A02
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("AMUC", "Round4", tabFlagAmucR4)

AMUC, Round 3

tabFlagAmucR3 <-
"
Strain Round plateType WellID
AMUC Round3 wBCAA A01
AMUC Round3 wBCAA A02
AMUC Round3 wBCAA A03
AMUC Round3 wBCAA B01
AMUC Round3 wBCAA B02
AMUC Round3 wBCAA B03
AMUC Round3 wBCAA A04
AMUC Round3 wBCAA A05
AMUC Round3 wBCAA B06
AMUC Round3 Plain A03
AMUC Round3 Plain B03
AMUC Round3 Plain A09
AMUC Round3 Plain G03
AMUC Round3 Plain G04
AMUC Round3 Plain G05

" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("AMUC", "Round3", tabFlagAmucR3)

AMUC, Round 2. The impacted growth issues in Round 2 positive controls give pause, but helpful to show here, nevertheless.

tabFlagAmucR2 <-
"
Strain Round plateType WellID
AMUC Round2 Plain E01
AMUC Round2 Plain A02
AMUC Round2 Plain B02
AMUC Round2 Plain E04
AMUC Round2 Plain D08
AMUC Round2 Plain B08
AMUC Round2 Plain E08
AMUC Round2 Plain A07
AMUC Round2 Plain C07
AMUC Round2 Plain E07
AMUC Round2 Plain B07
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("AMUC", "Round2", tabFlagAmucR2)

CBUT curation by round

CBUT, Round 4

tabFlagCbutR4 <-
"
Strain Round plateType WellID
CBUT Round4 wBCAA D12
CBUT Round4 Plain C01
CBUT Round4 Plain B01
CBUT Round4 Plain B03
CBUT Round4 Plain B04
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("CBUT", "Round4", tabFlagCbutR4)

CBUT, Round 3. Positive controls weirdly lower yield than usual. This is a challenging run for CBUT. Plucking out the most egregious anomalies, anyway. Looks like the impact growth + SFU was a pretty severe insult. BCAA didn’t rescue anything here, either.

tabFlagCbutR3 <-
"
Strain Round plateType WellID
CBUT Round3 Plain A02
CBUT Round3 Plain A03
CBUT Round3 Plain A07
CBUT Round3 Plain A09
CBUT Round3 Plain E04
CBUT Round3 Plain C07
CBUT Round3 wBCAA A05
CBUT Round3 wBCAA A07
CBUT Round3 wBCAA A08
CBUT Round3 wBCAA B01
CBUT Round3 wBCAA B02
CBUT Round3 wBCAA B03
CBUT Round3 wBCAA B04
CBUT Round3 wBCAA B05
CBUT Round3 wBCAA B06
CBUT Round3 wBCAA B07
CBUT Round3 wBCAA B08
CBUT Round3 wBCAA B09
CBUT Round3 wBCAA C07
CBUT Round3 wBCAA C08
CBUT Round3 wBCAA C09
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("CBUT", "Round3", tabFlagCbutR3)

CBUT, Round 2. This looks respectably clean. The effect appears to be mostly noticeable at the highest concentration in this smaller titration. Looks like a couple controls had an issue, but otherwise fine.

# CBUT Round2 Plain G07
# CBUT Round2 Plain F09
# CBUT Round2 Plain F10
tabFlagCbutR2 <-
"
Strain Round plateType WellID
CBUT Round2 Plain G09
CBUT Round2 Plain F11
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("CBUT", "Round2", tabFlagCbutR2)

CBUT, Round 1. Some early-clipped curves could affect the lag calculation, but the aggregated curve would be fine and the sd between replicates look good. This spans a range of nearly no effect to clearly pretty big effect, despite the relatively large impact of 3% DMSO.

tabFlagCbutR1 <-
"
Strain Round plateType WellID
CBUT Round1 Plain B03
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("CBUT", "Round1", tabFlagCbutR1)

CBEI curation by rounds

CBEI, Round 4

tabFlagCbeiR4 <-
"
Strain Round plateType WellID
CBEI Round4 Plain B04
CBEI Round4 Plain G04
CBEI Round4 Plain C01
CBEI Round4 wBCAA F09
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("CBEI", "Round4", tabFlagCbeiR4)

CBEI, Round 3

tabFlagCbeiR3 <-
"
Strain Round plateType WellID
CBEI Round3 Plain G07
CBEI Round3 Plain F01
CBEI Round3 Plain A01
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("CBEI", "Round3", tabFlagCbeiR3)

CBEI, Round 2

tabFlagCbeiR2 <-
"
Strain Round plateType WellID
CBEI Round2 Plain A07
CBEI Round2 Plain B07
CBEI Round2 Plain E07
CBEI Round2 Plain B02
CBEI Round2 Plain E02
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("CBEI", "Round2", tabFlagCbeiR2)

EHAL curation by round

EHAL, Round 5

tabFlagEhalR5 <-
"
Strain Round plateType WellID
EHAL Round5 Plain E10
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("EHAL", "Round5", tabFlagEhalR5)

plotly::ggplotly()

BINF curation by round

BINF, Round 5

tabFlagBinfR5 <-
"
Strain Round plateType WellID
BINF Round5 Plain A01
BINF Round5 Plain A04
BINF Round5 Plain A10
BINF Round5 Plain A12
" %>% 
  fread(stringsAsFactors = FALSE)
plot_flag("BINF", "Round5", tabFlagBinfR5)

plotly::ggplotly()

Whale plots

Define whale plot function, to simplify and enforce consistency plot summaries.

plot_whale = function(tab1, 
                      title = "", 
                      rangeOrSmooth = "range", 
                      smoothSpan = 0.5, 
                      seLevel = 0.99,
                      ctlMaxHrsDmso = Inf,
                      ctlMaxHrsWater = Inf,
                      alphaWhale = 0.25){
  p1 <- tab1WaterOnly <- tab1DmsoOnly <- tab1ConcLabel <- NULL
  tab1WaterOnly <-
    tab1 %>% copy %>% 
    # Omit the water-only (no DMSO) control from main data
    .[(Solvent == "Water")] %>% 
    # Only control
    .[(Solute == "None")] %>% 
    # Set to NULL so that it shows up everywhere
    .[, Solute := NULL] %>% 
    .[, Titration := NULL] %>% 
    # Trim off weirdness later on in growth curve
    .[(Hours <= ctlMaxHrsWater)] %>%
    # Coarse-grain hours to avoid some weird effects.
    .[, Hours := round(Hours)] %>% 
    # Compute median, min, max at each hour
    .[, .(odMed = median(OdSplineBaselineSubtract), 
          odMin = min(OdSplineBaselineSubtract), 
          odMax = max(OdSplineBaselineSubtract)), 
      by = .(Strain, Round, plateType, Hours)]
  tab1DmsoOnly <-
    tab1 %>% copy %>% 
    # Omit the water-only (no DMSO) control from main data
    .[(Solvent == "DMSO")] %>% 
    # Only control
    .[(Solute == "None")] %>% 
    # Set to NULL so that it shows up everywhere
    .[, Solute := NULL] %>% 
    .[, Titration := NULL] %>% 
    # Trim off weirdness later on in growth curve
    .[(Hours <= ctlMaxHrsDmso)] %>%
    # Coarse-grain hours to avoid some weird effects.
    .[, Hours := round(Hours)] %>%
    # Compute median, min, max at each hour
    .[, .(odMed = median(OdSplineBaselineSubtract), 
          odMin = min(OdSplineBaselineSubtract), 
          odMax = max(OdSplineBaselineSubtract)), 
      by = .(Strain, Round, plateType, Hours)]

  # Table for inscribing the SFU concentration on each facet
  tab1ConcLabel <-
    tab1 %>% copy %>% 
    # Filter controls to max hours
    .[!(Solute == "None" & Hours > ctlMaxHrsWater)] %>% 
    # Set dummy Hours and OD values for positioning label
    .[, Hours := 0.0] %>% 
    .[, OdSplineBaselineSubtract := max(OdSplineBaselineSubtract)] %>% 
    # Now that max calculated, filter out controls completely
    .[(Solute != "None")] %>% 
    # downselect to take unique-entry table
    .[, .(Strain, Round, SolventPercent, plateType,
          Solute, Titration, ConcentrationMillimolar,
          Hours, OdSplineBaselineSubtract)] %>% 
    unique() %>% 
    # text for printing on panel
    # Solute Concentration
    .[, showConc := round(ConcentrationMillimolar, digits = 2)] # %>% paste0(., " [mM]\nDMSO: ", SolventPercent, "%")]
  
  # Define plot
  p1 <-
    tab1 %>%
    # Omit the controls for now (will highlight region on all facets)
    .[(Solute != "None")] %>% 
    setorder(Strain, Round, plateType, WellID, Hours) %>% 
    ggplot(aes(Hours, OdSplineBaselineSubtract, 
               color = plateType)) + 
    facet_grid(
      rows = vars(Titration),
      cols = vars(Solute)
    )
  ## Add control layers
  # DMSO control
  if(rangeOrSmooth == "range"){
    p1 <- p1 +
      geom_ribbon(
        data = tab1DmsoOnly,
        color = NA,
        fill = "darkgreen",
        size = 0.1,
        alpha = alphaWhale,
        mapping = aes(y = odMed, ymin = odMin, ymax = odMax)
      )
  } else {
    p1 <- p1 +
      stat_smooth(
        geom = "ribbon",
        color = NA,
        span = smoothSpan,
        level = seLevel,
        data = tab1 %>% copy %>% 
          # Omit the water-only (no DMSO) control from main data
          .[(Solvent == "DMSO")] %>% 
          # Only control
          .[(Solute == "None")] %>% 
          # Trim off weirdness later on in growth curve
          .[(Hours <= ctlMaxHrsDmso)] %>%
          # Set to NULL so that it shows up everywhere
          .[, Solute := NULL] %>% 
          .[, Titration := NULL],
        se = TRUE,
        fill = "darkgreen",
        alpha = alphaWhale
      )
  }
  # No-DMSO / 'Water' control
  if(rangeOrSmooth == "range"){
    p1 <- p1 +
      geom_ribbon(
        data = tab1WaterOnly,
        color = NA,
        fill = "darkblue",
        size = 0.1,
        alpha = alphaWhale,
        mapping = aes(y = odMed, ymin = odMin, ymax = odMax)
      )
  } else {
    p1 <- p1 +
      stat_smooth(
        geom = "ribbon",
        color = NA,
        span = smoothSpan,
        level = seLevel,
        data = tab1 %>% copy %>% 
          # Omit the water-only (no DMSO) control from main data
          .[(Solvent == "Water")] %>% 
          # Only control
          .[(Solute == "None")] %>% 
          # Trim off weirdness later on in growth curve
          .[(Hours <= ctlMaxHrsWater)] %>%
          # Set to NULL so that it shows up everywhere
          .[, Solute := NULL] %>% 
          .[, Titration := NULL],
        se = TRUE,
        fill = "darkblue",
        alpha = alphaWhale
      )
  }
  p1 <- p1 +
    # "the data"
    geom_path(
      mapping = aes(
        group = paste(Strain, Round, plateType, WellID), 
        y = OdSplineBaselineSubtract),
      alpha = 1, 
      size = 0.5) +
    scale_color_manual(values = c(Plain = "black", wBCAA = "brown")) +
    # Add label for concentration of this panel row
    geom_text(
      data = tab1ConcLabel,
      color = "black",
      mapping = aes(
        y = OdSplineBaselineSubtract,
        label = showConc),
      hjust = 0,
      vjust = 1,
      size = 2.2,
      alpha = 0.5
    ) +
    # scale_color_brewer(palette = "Set2") +
    ylab("OD (spline-smooth, baseline subtracted)") +
    # Set a consistent number of y-axis breaks
    scale_y_continuous(
      breaks = scales::breaks_pretty(n = 2)
      # breaks = scales::breaks_extended(2, m = 2)
      # breaks = scales::breaks_extended(2)
    ) + 
    theme(
      legend.position = "none",
      strip.text.y = element_blank()
    ) +
    theme(panel.grid = element_blank())
  if(title != ""){
    p1 <- p1 + ggtitle(label = title)
  }
  return(p1)
}

filter function

filter_whale = function(tabGm, iStrain, iRound){
  tabWhale <-
    tabGm %>% copy %>% 
    .[(Inoculation == "Inoculated")] %>% 
    .[(Strain == iStrain)] %>% 
    .[(Solute != c("Sulfadiazine"))] %>% 
    .[(Round == iRound)] %>% 
    .[, ConcRound := round(Concentration, digits = 2)] %>% 
    setorder(Strain, Round, plateType, WellID, Hours)
  return(tabWhale)
}

Whale plots, AMUC

plot_whale(
  tabGrowthMask %>% 
    filter_whale("AMUC", "Round2") %>% 
    .[!tabFlagAmucR2,
      on = c("Strain", "Round", "plateType", "WellID")], 
  title = "AMUC, Round 2")

pWhaleAmucR3 <- 
  plot_whale(
    tabGrowthMask %>% 
      filter_whale("AMUC", "Round3") %>% 
      .[!tabFlagAmucR3, 
        on = c("Strain", "Round", "plateType", "WellID")] %>% 
      .[(plateType == "Plain")],
    rangeOrSmooth = "range",
    # smoothSpan = 3, seLevel = 0.99999,
    ctlMaxHrsDmso = 50, 
    ctlMaxHrsWater = 38,
    title = "AMUC, Round 3")
pWhaleAmucR3

plot_whale(
  tabGrowthMask %>% 
    filter_whale("AMUC", "Round4") %>% 
    .[!tabFlagAmucR4,
      on = c("Strain", "Round", "plateType", "WellID")] %>%
    .[(plateType == "Plain")],
  ctlMaxHrsDmso = 45,
  ctlMaxHrsWater = 38, 
  title = "AMUC, Round 4")

Whale plots, CBUT

plot_whale(
  tabGrowthMask %>% 
    filter_whale("CBUT", "Round1") %>% 
    .[!tabFlagCbutR1,
      on = c("Strain", "Round", "plateType", "WellID")],
  title = "CBUT, Round 1")

plot_whale(
  tabGrowthMask %>% 
    filter_whale("CBUT", "Round2") %>% 
    .[!tabFlagCbutR2,
      on = c("Strain", "Round", "plateType", "WellID")], 
  title = "CBUT, Round 2",
  ctlMaxHrsDmso = 15,
  ctlMaxHrsWater = 15,
  rangeOrSmooth = "smooth",
  smoothSpan = 0.5)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

plot_whale(
  tabGrowthMask %>% 
    filter_whale("CBUT", "Round3") %>%
    .[!tabFlagCbutR3,
      on = c("Strain", "Round", "plateType", "WellID")] %>%
    # drop BCAA condition for plotting purposes
    .[(plateType == "Plain")],
  ctlMaxHrsDmso = 25,
  ctlMaxHrsWater = 20,
  title = "CBUT, Round 3",
  rangeOrSmooth = "smooth",
  smoothSpan = 0.5)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

plot_whale(
    tabGrowthMask %>% 
    filter_whale("CBUT", "Round4") %>%
    # Counter-select the flagged wells
    .[!tabFlagCbutR4,
      on = c("Strain", "Round", "plateType", "WellID")] %>%
    # drop BCAA condition for plotting purposes
    .[(plateType == "Plain")],
  ctlMaxHrsDmso = 25,
  ctlMaxHrsWater = 20,
  title = "CBUT, Round 4",
  rangeOrSmooth = "smooth",
  smoothSpan = 0.5)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Whale plots, CBEI

plot_whale(
  tabGrowthMask %>% 
    filter_whale("CBEI", "Round2") %>% 
    .[!tabFlagCbeiR2,
      on = c("Strain", "Round", "plateType", "WellID")],
  title = "CBEI, Round 2")

plot_whale(
  tabGrowthMask %>% 
    filter_whale("CBEI", "Round3") %>% 
    # drop BCAA condition for plotting purposes
    .[(plateType == "Plain")] %>% 
    .[!tabFlagCbeiR3,
      on = c("Strain", "Round", "plateType", "WellID")],
  ctlMaxHrsDmso = 19,
  title = "CBEI, Round 3")

plot_whale(
  tabGrowthMask %>% 
    filter_whale("CBEI", "Round4") %>% 
    # drop BCAA condition for plotting purposes
    .[(plateType == "Plain")] %>% 
    .[!tabFlagCbeiR4,
      on = c("Strain", "Round", "plateType", "WellID")],
  title = "CBEI, Round 4")

Whale plot, EHAL

plot_whale(
  tabGrowthMask %>% 
    filter_whale("EHAL", "Round5") %>% 
    # No BCAA condition, showing this explicitly
    .[(plateType == "Plain")] %>% 
    .[!tabFlagEhalR5,
      on = c("Strain", "Round", "plateType", "WellID")],
  title = "EHAL, Round 5")

Whale plot, BINF

plot_whale(
  tabGrowthMask %>% 
    filter_whale("BINF", "Round5") %>% 
    # No BCAA condition, showing this explicitly
    .[(plateType == "Plain")] %>% 
    .[!tabFlagBinfR5,
      on = c("Strain", "Round", "plateType", "WellID")],
  title = "BINF, Round 5")

Pick representative slices for main figure

This entire report is being shared as supplementary figure, for reference and to document the data censoring choices. The goal of slicing here is to pick a representative combination that accurately communicates the inhibition (or lack thereof) without extraneous idiosyncracies that distract from the relative ranking and qualitative outcome.

Also, only so many panels will fit in a main figure plot.

pWhaleMainAmuc <-
  plot_whale(
    tabGrowthMask %>% 
      filter_whale("AMUC", "Round3") %>% 
      .[!tabFlagAmucR3, 
        on = c("Strain", "Round", "plateType", "WellID")] %>% 
      .[(plateType == "Plain")] %>% 
      .[(Solute == "Glimepiride" & Titration %in% c(2, 3, 4) |
           Solute == "Glipizide" & Titration %in% c(2, 3, 4) |
           Solute == "Glyburide" & Titration %in% c(2, 3, 4) |
           Solute == "None")], 
    rangeOrSmooth = "range",
    ctlMaxHrsDmso = 50, 
    ctlMaxHrsWater = 38,
    title = "AMUC")
pWhaleMainAmuc

pWhaleMainCbut <-
  plot_whale(
    tabGrowthMask %>% 
      filter_whale("CBUT", "Round1") %>% 
      .[!tabFlagCbutR1, 
        on = c("Strain", "Round", "plateType", "WellID")] %>% 
      .[(plateType == "Plain")] %>% 
      .[(Solute == "Glimepiride" & Titration %in% c(2, 3, 4) |
           Solute == "Glipizide" & Titration %in% c(2, 3, 4) |
           Solute == "Glyburide" & Titration %in% c(1, 2, 3) |
           Solute == "None")] %>% 
      # update titration rank 
      # so that there are not extra dummy rows in the final plot
      .[(Solute == "Glyburide"), Titration := Titration + 1L],
    rangeOrSmooth = "range",
    title = "CBUT")
pWhaleMainCbut

pWhaleMainCbei <-
  plot_whale(
    tabGrowthMask %>% 
      filter_whale("CBEI", "Round3") %>% 
      .[!tabFlagCbutR1, 
        on = c("Strain", "Round", "plateType", "WellID")] %>% 
      .[(plateType == "Plain")] %>% 
      .[(Solute == "Glimepiride" & Titration %in% c(1, 2, 3) |
           Solute == "Glipizide" & Titration %in% c(1, 2, 3) |
           Solute == "Glyburide" & Titration %in% c(1, 2, 3) |
           Solute == "None")], 
      # .[, Titration := 0],
    rangeOrSmooth = "range",
    ctlMaxHrsDmso = 19, 
    ctlMaxHrsWater = Inf,
    title = "CBEI")
pWhaleMainCbei

pWhaleMainBinf <-
  plot_whale(
    tabGrowthMask %>% 
      filter_whale("BINF", "Round5") %>% 
      .[!tabFlagBinfR5, 
        on = c("Strain", "Round", "plateType", "WellID")] %>% 
      .[(Hours < 18)] %>% 
      .[(plateType == "Plain")] %>% 
      .[(Solute == "Glimepiride" & Titration %in% c(3, 4, 5) |
           Solute == "Glipizide" & Titration %in% c(3, 4, 5) |
           Solute == "Glyburide" & Titration %in% c(3, 4, 5) |
           Solute == "None")], 
      # .[, Titration := 0],
    rangeOrSmooth = "range",
    ctlMaxHrsDmso = 17, 
    ctlMaxHrsWater = 17,
    title = "BINF")
pWhaleMainBinf

pWhaleMainEhal <-
  plot_whale(
    tabGrowthMask %>% 
      filter_whale("EHAL", "Round5") %>% 
      .[!tabFlagEhalR5, 
        on = c("Strain", "Round", "plateType", "WellID")] %>% 
      .[(plateType == "Plain")] %>% 
      .[(Solute == "Glimepiride" & Titration %in% c(3, 4, 5) |
           Solute == "Glipizide" & Titration %in% c(3, 4, 5) |
           Solute == "Glyburide" & Titration %in% c(3, 4, 5) |
           Solute == "None")], 
      # .[, Titration := 0],
    rangeOrSmooth = "range",
    ctlMaxHrsDmso = Inf, 
    ctlMaxHrsWater = Inf,
    title = "EHAL")
pWhaleMainEhal

Build Figure 5

layout <- "
AAABBBBB
AAABBBBB
AAABBBBB
CCCCCCCC
CCCCCCCC
CCCCCCCC
"
whalePlotStripSize <- 6
whalePlotTitleSize <- 6
whalePlotAxisSize <- 6
pFig5 <- NULL
pFig5 <-
  (
    pSfuA1c +
      theme(
        # plot.title = element_text(hjust = 0, size = 5),
        strip.text = element_text(hjust = 0, size = 6),
        axis.title = element_blank(),
        axis.text.x = element_text(hjust = 0.5, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6)
      )
  ) + 
  (
    pTestsBwGrp + 
      theme(
        plot.title = element_blank(), 
        axis.title = element_blank(),
        axis.text.x = element_text(hjust = 0.5, size = 6),
        axis.text.y = element_text(hjust = 1, size = 6),
        strip.text = element_text(hjust = 0, size = 6)
      )
  ) +
  (
    (pWhaleMainAmuc + 
       ylab("OD\n600 nm") +
       theme(
         title = element_text(size = whalePlotTitleSize),
         strip.text = element_text(size = whalePlotStripSize),
         axis.text = element_text(size = whalePlotAxisSize),
         axis.title.y = element_text(
           margin = margin(r = 10, unit = "pt"),
           angle = 0, 
           hjust = 1, 
           vjust = 1, 
           size = 6), 
         axis.title.x = element_blank()
       )
    ) +
      (pWhaleMainBinf + 
         theme(
           title = element_text(size = whalePlotTitleSize),
           strip.text = element_text(size = whalePlotStripSize),
           axis.text = element_text(size = whalePlotAxisSize),
           axis.title.y = element_blank(), 
           axis.title.x = element_blank())
      ) +
      (pWhaleMainCbei + 
         theme(
           title = element_text(size = whalePlotTitleSize),
           strip.text = element_text(size = whalePlotStripSize),
           axis.text = element_text(size = whalePlotAxisSize),
           axis.title.y = element_blank(), 
           axis.title.x = element_blank()
         )
      ) +
      (pWhaleMainCbut + 
         theme(
           title = element_text(size = whalePlotTitleSize),
           strip.text = element_text(size = whalePlotStripSize),
           axis.text = element_text(size = whalePlotAxisSize),
           axis.title.y = element_blank(),
           axis.title.x = element_blank())
      ) +
      (pWhaleMainEhal + 
         theme(
           title = element_text(size = whalePlotTitleSize),
           strip.text = element_text(size = whalePlotStripSize),
           axis.text = element_text(size = whalePlotAxisSize),
           axis.title.y = element_blank(), 
           axis.title.x = element_text(size = 6)
         )
      ) +
      plot_layout(nrow = 2)
  ) +
  plot_layout(design = layout) +
  plot_annotation(tag_levels = list(letters[1:3])) &
  theme(plot.tag = element_text(size = 8, face = "bold", family = "Sans"))
# pFig5
# Journal guided formatting, dimensions, font sizes
figHeight = 225

Save to standalone figure files.

ggsave("Figure-05.pdf", pFig5, 
       device = cairo_pdf,
       width = 170,
       # Have room up to 200 or so?
       height = figHeight, 
       dpi = 300, 
       units = "mm")